%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%% Solve for Optimal Contract %%%
%Applies in general case with maturity delta>0
%Baseline is obtained via delta=0
%The procedure is:
% 1) Solve the model for given q
% 2) Iterate over q
% 3) Choose the optimal q that leads to highest surplus
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

qSpace=linspace(0,qBar,qGrid);

for z=1:qGrid
    
    q=qSpace(z);
    
    %Solve in SB.m the model without screening moral hazard, i.e., the
    %second best case, to obtain solution at boundary V^B
 
    SB;
    %Boundary conditions for F' at V^B
    dFStart=0;
    
    % Slope at the boundary V^B is zero if and only if aB
    % For aB=0, different boundary condition applies

    if (aB<0.00000001)

        dFStart=(FB-phi*(gamma-r))/phi;
        
        
    end
    
    %right indicates whether V_0=kappa*q>VB

    right=1;
    
    if (kappa*q<VB)
        
        right=0;
        
    end
    
    %Solve the HJB equation starting from V^B to V_0
    %F'(V)=0 in the limit V\to V^B, we intialize the slope with a small
    %negative value if V^0>V^B and a^B>0
    %We pick tolerance=0.0001; see para.m file

    sol = ode15s(@ffunODE,[VB max(kappa*q,VB+0.001)+0.001*(kappa*q==VB)],[FB dFStart+tolerance-2*tolerance*right],options, r, gamma ,Lambda, phi,q,aBar,delta);

    
    %Get starting value F_0
    F=sol.y(1,end);
    dF=sol.y(2,end);
    aEnd=(F-dF*(sol.x(end)+phi)-phi*(gamma-r))/phi;
    WBound=(sol.y(1,1)-dFStart*(sol.x(1)+phi)-(gamma-r)*phi);
    
    
    
    [m,ind]=min(abs(sol.x-kappa*q));
    V0=sol.x(ind);
    F0=sol.y(1,ind);
    
    %Collect outputs in vectors
    %Calculate F_{0^-}
    vF0(z)=F0-0.5*kappa*q^2;
    vAB(z)=aB;
    vAEnd(z)=aEnd;
    vV0(z)=V0;
    
    %Break the loop when q leads to degenerate outcome
    %In this case, we stop the iteration, as higher q then also would lead
    %to degenerate outcomes
    if (vF0(z)<-10)

        vF0(z+1:qGrid)=zeros(qGrid-z,1);

        break;

    end
    
    
end

%Determine Optimal q
[m,ind]=max(vF0);
q=qSpace(ind);
a0=vAEnd(ind);
V0=vV0(ind);

